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ABSTRACT 

A theoretical investigation o£ the effects of short period 
angle of attack oscillation on air-to-ground weapons dive 
delivery accuracy was conducted. The investigation included 
the computer modeling of the lonjitudinal control system of 
the F-4 aircraft. This longitudinal control system was used 
to fly the aircraft through its standard conventional dive 
bomb delivery maneuver, making corrections for various dive 
angle errors, detected four thousand feet above release 
altitude . 

The dynamics of the F-4 aircraft motion was represented 
by a matrix of dimensional, body axis longitudinal stability 
derivatives and force coefficients. This matrix representa- 
tion of the aircraft dynamics was used to solve for the time 
history of aircraft perturbation, velocity, angle of attack, 
pitch angle, pitch rate and flight path angle, following 
correction to the desired delivery maneuver flight path. 

The perturbation angle of attack is a measure of the 
gunsight aimpoint error at any given time, if all other de- 
livery parameters are met, and converts directly to weapon 
impact error. 
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DEFINITIONS OF SYMBOLS 



REFERENCE QUANTITIES 

Y 

To - total aircraft velocity in reference state 

Uq - aircraft velocity parallel to FRL in 

reference state 

Wq - aircraft velocity perpendicular to FRL 

in reference state 

Mn - aircraft reference Mach number 
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ith 

FRL 



reference flight path angle 
reference pitch angle 
reference angle of attack 
aircraft mass 
aircraft weight 
air density 
wing area 

wing mean aerodynamics chord 

dynamic pressure 

acceleration of gravity 

moment of inertia about y axis 

thrust incidence with FRL + 

perpendicular distance to thrust line 
from c.g. 

longitudinal reference axis -fuselage 
reference line 

longitudinal distance from c.g. to 
bob weight 



(ft/sec) 

(ft/sec) 

(ft/sec) 



(radians) 
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(pounds) 

(slugs/ft^) 
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PERTURBATION QUANTITIES 



Y - perturbation flight path angle (rad or deg) 

0 - perturbation pitch angle (rad) 

a - perturbation angle of attack (rad/milli- 

radians 
when refer 
to sight 
angle) 

- perturbation velocity component along FRL (ft/sec) 
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q 
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- perturbation velocity component perpen- 
dicular to FRL 

- perturbation pitch rate 

- perturbation stabilator position 



TOTAL QUANTITIES 



0 

ALPHA 

U 

M 



- total flight path angle 
r= Tj + Y = 0 - ALPHA 

- total pitch angle 0 = 0. + 0 



- total FRL angle of attack 
ALPHA = ttj + a 

- total velocity compert along FRL 
U = U„ + u 

- total Mach number 



(ft/ sec) 

(radians/sec) 

(rad-pos. 

L.E. up) 



(rad/ deg) 
(rad) 

(rad/mils) 
(ft/ sec) 



STABILITY AXIS DIMENSIONLESS COEFFICIENTS AND DERIVATIVES 
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All moment coefficients and derivatives are the same for body- 
axis as for stability axis. 
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I. INTRODUCTION 



The theory and practice of fixed gunsight, air-to-ground 
weaponry is well established and differs only slightly from 
one aircraft to another due to minor handling quality and 
weapons release system differences among aircraft. Just as 

I 

well established is the analysis of impact error sensitivi- 
ties to variation of actual from desired release parameters. 
Table I shows the desired release parameters and error sen- 
sitivities in the most widely used F-4 dive bomb delivery. 

It is not the purpose of this study to disagree with any of 
the release parameters for this or any other delivery maneu- 
ver, but rather to demonstrate that increased emphasis in 
training on one of the error sensitivity categories, that of 
gunsight placement as influenced by FRL angle of attack per- 
turbations, can greatly increase pilot competence in air-to- 
ground weapons delivery. 

A lack of emphasis on angle of attack input to gunsight 
picture seems justified on the surface, since everything 
which is an input to sight angle can be computed and set 
prior to the delivery maneuver and should not be subject to 
error. In addition, Table I shows that the impact error 
due to a sight picture error on one mil (0.1% of a radian) 
is smaller than the impact error due to an airspeed error 
of 10 knots, or a flight path dive angle error of one de- 
gree. For these reasons and since only one FRL angle of 
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Table I. 



Desired Delivery Parameters and Error 
Sensitivities for Maneuver Studied 



Release Parameters 



Dive Angle 


40 Degrees 


Release Airspeed 


500 KTAS 


Release Height 


5000 Feet AGL 


Sight Angle 


•116 Mils 


Release Horiz Range 


4720 Feet 


Release Slant Range 


6875 Feet 


Parameter Error 


Impact Error 


Dive Angle 


1 Steep 


33 Feet long 


1° Shallow 


35 Feet short 


Release Height 


100 Feet High 


13 Feet short 


100 Feet Low 


13 Feet long 


Release Airspeed 


10 KTS Fast 


35 Feet long 


10 KTS Slow 


37 Feet short 


Sight Angle/Aim Point 


+ 1 Mil/1 Mil Long 


10 Feet long 


- 1 Mil/1 Mil Short 


10 Feet short 



(This information was calculated from data in Ref. 3. Sight 



angle includes weapon trajectory drop, FRL 



angle of attack 



and gunsight parallax correction.) 



attack is associated with one steady flight condition at a 
given weight and configuration, it has been accepted that the 
angle of attack input to the sight setting is a quantity which 
varies only with release airspeed and dive angle. This is 
true, however, only if the aircraft is free from all longi- 
tudinal oscillations rising from the most recent atmospheric 
disturbance or if longitudinal control action were taken. 
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Most of the emphasis given to sight error due to angle of 
attack perturbations in the past can be summed up in the 
statement: "The aircraft must be in one-G flight at release 

in order for pipper placement to be an accurate indication 
of bomb impact point." Very little emphasis has been placed 
on the fact that a finite time interval must expire after 
the most recent longitudinal control input before the actual 
FRL angle of attack matches the steady state angle of attack 
for the desired flight condition. If bomb release occurs 
during this finite interval, any difference between perturbed 
FRL angle of attack and angle of attack computed for the 
sight setting will manifest itself as a bomb impact error. 
Finding the magnitude and duration of this angle of attack 
perturbation after making various dive angle error correc- 
tions was the purpose of this investigation. 
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II. PROBLEM SET-UP 



The standard F-4 visual dive delivery of the MK 82 
LDGP Bomb with conical fins was chosen for the investiga- 
tion. This maneuver involves a roll-in from 12,000 feet at 
400 KTAS on a heading 90° to desired run-in line. Roll-in 
is accomplished by a 135° bank angle descending turn to a 
40° dive along the run-in line. Roll-in is accomplished at 
military power and completed by approximately 10,000 feet. 
Power is reduced to idle at 430 KTAS (at approximately 9,000 
feet) , and the aircraft accelerates at idle power in the 40° 
dive to 500 KTAS at the release altitude of 5,000 feet AGL. 

The largest impact error in range is due to dive angle 
error at release. The causes for this error, shown in 
Figure 1, are: 

1. In a dive angle other than 40°, with the gunsight 
pipper on the target, the range to the target is different 
than that planned for at 5,000 feet, while the bomb range 
from this point is approximately the same as if in the plan- 
ned dive angle . 

2. The FRL angle of attack is dependent on dive angle, 
but the input to the gunsight for angle of attack is based 
on a 40° dive angle. 

In this study, pilot detection of a dive angle error at 
9,000 feet (a typical point to detect an error, shortly 
after roll-out on final run in line) caused the pilot to 
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ww 





make an additional dive angle change to correct for the error 
in dive angle. 

There are two different situations which can be considered 
to be the same dive angle error. A five degree shallow error 
exists for either of the following cases: 

1. The pipper placement is correct for 9,000 feet alti- 
tude but the aircraft dive angle, as noted visually outside 
the aircraft or on the aircraft attitude indicator, is 35°. 

2. The aircraft is actually in a 40° dive but the pipper 
placement at 9,000 feet altitude is 5° short of where it 
should be on the ground at that point. 

In either case the aircraft is approximately 1,850 feet 
(horizontal range) farther from the target than it should be 
at that altitude in the maneuver, and the pilot must shallow 
the dive still further to re-position the aircraft on the 
proper flight path over the ground, and then re-establish 
the desired 40° dive angle. Since pilots rapidly acquire 
skill at rolling into a desired dive angle, the normal dive 
angle error situation is to be in a 40° dive but find the 
pipper aimed too far short or not short enough of the target 
at the 9,000 feet altitude check point. This was the case 
considered in this investigation. 

A. DIFFERENCES BETWEEN ACTUAL MANEUVER AND COMPUTER SIMULATION 

1 . Actual Maneuver 

In an actual run the aircraft is constantly increasing 
airspeed from approximately 430 KTAS at 9,000 feet to 500 KTAS 
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at 5,000 feet with changing lift and drag coefficients. Air 
density increases as altitude decreases and typically target 
elevation is not at sea level. 



2 . Computer Simulation 

In the computer simulation it was assumed that re- 
lease airspeed was attained at 9,000 feet and no further 
acceleration was achieved during the runs except as a result 
of perturbations due to longitudinal maneuvering. Since the 
assumed quasi-steady state velocity is the actual desired 
release velocity, results in the area of release altitude 
were very accurate. Air density was assumed constant at 
standard day density for 7,000 feet M.S.L., the mid-altitude 
for the run. Target elevation was assumed to be standard 
sea level. Other assumptions and conventions used were: 

a. Flat earth fixed in space, 

b. Atmosphere fixed to earth (no wind), 

c. Rigid airframe, 

d. Body axis system along FRL used throughout, 

e. X-Z plane of symmetry, and 

f. Small perturbation angles and velocities. 



B. REPRESENTATION OF AIRFRAME DYNAMICS 



The following body axis force equation system was ob- 
tained from Ref. 1 in the Laplace domain. 
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u u 
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a a q 




0 

■> — ‘ 




. sj 



19 



where notation and definitions are given in Appendix A. The 
above system is represented in the time domain by the follow- 
ing linear differential equations. 



1 . 



2 . 



u = (X^ + cosOu + X^a - gcosG^G - (W, - X^)g + X^ 65 
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4 . q = (M + 7 T)u + M*a + Ma + Mq + Mj(. 6 
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Substitution of Equation 2 for a into equation 4 yields 
the following time domain matrix differential equation in 
the perturbation quantities u, a, 0 and q. 
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C. COMPUTER REPRESENTATION OF F-4 LONGITUDINAL CONTROL SYSTEM 



Reference 1 contains diagrammatic representation of the 
individual components of the F-4 longitudinal control system 
as a linear control feedback system which can be combined to 
form the complete linear control system shown in Figure 2. 
Either Block Diagram Algebra or the Signal Flow Graph tech- 
nique may be used to reduce the system in Figure 2 to the 
simplified closed loop system shown in matrix equation form 
in Figure 3. Derivation by Block Diagram Algebra techniques 
of the following Laplace domain representation of the con- 
trol system is presented in detail in Appendix A. 

(l+G2G^G^HlH2(S) - G4G5H3H,(S))6g(S) = G i G^G sG^ (S) [ - Y (S) ] . 

As shown in Appendix A, precise representations in the 
Laplace domain of all the control system transfer functions 
would lead to a 42nd order differential equation, the numer- 
ical solution to which would be impractical for the purposes 
of this study. For this reason a method was sought to sim- 
plify the more complicated transfer functions which were not 
directly indicative of airframe dynamics: the aircraft feel 

system, defined GaCS); the stabilator actuator, defined Gi» (S) ; 
and the vertical acceleration feedback system to the Bob- 
weight, defined Hi(S). These simplifications are shown in 
Appendix A. The 0/6g(S) transfer function, defined Gs(S), 
although complicated, was not simplified since it represented 
actual airframe dynamic response to stabilator movement and 
its accuracy was considered critical to the investigation of 
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Figure 2. F-4 Longitudinal Control System. 
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the problem. With the above simplifications, the following 
equation is the Laplace domain representation for the sta- 
bilator position in terms of pilot stick force and dive 
angle change required (y^"Y). ♦ (See Appendix A for deriva- 
tion.) 

(S® + 18.45S” + 35.57S^ + 95. OS^ + 1.201S + 0.0224) 6g(S) = 

- 0.01113Kp(S^ + 3.8025“ + 17.995^ + 16.375^ + 1.201S 

+ 0.0224) [y^-y(S)] . 

In the time domain this becomes the following fifth 
order linear differential equation: 

+ 18.456g^^^ + 35.576g^^^ + 95. 05^ + 1.2016g + 0.02246g 
- 0.01113Kp[(Y^-Y) + 3.802 (y^-Y) + 17 . 99 (y^- Y) 

+ 16.37(y^-y) + 1.201(y^-y) + 0 . 0224 (y^- Y) ] • 

The state space technique of Reference 5 for representa- 
tion of nth-order systems of linear differential equations, 
in which the forcing function involves derivative terms, was 
used to cast the control system as a set of five linear first 
order differential equations, in general, for the following 
fifth order linear differential equation, 

y^^^ + + a^y^^^ + a^y + agy + a^y = 

- + 3gU^^^ + Byii + BgU + BgU 

with the following definition of terms 
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z = {1 0 0 0 0} {x^x^x^XgXg + bg [u] . 

In the case of the control system, z is defined to be the 
stabilator position , and u is the forcing function (y "Y) . 

o C 

The subscripts 5 through 9 were used since the control sys- 
tem equations were to be coupled with the four airframe 
equations of motion to form a ninth-order system for numeri- 
cal solution. 

D. AIRFRAME AND CONTROL SYSTEM COUPLED 

The airframe and control system equation were coupled 

through defined by x_-K b„(y -y) . Using the fact that 

o D p U C 

flight path angle F is defined by F = 0-ALPHA and the flight 
path perturbation angle is defined by y = 9 - a, the final 
set of nine coupled linear differential equations appears as 
the feedback control system of Figure 3. 



26 



III. PROBLEM SOLUTION AND RESULTS 



A. CASES CONSIDERED 

Two separate initial conditions were treated. In the 
first case the pilot discovered at 9,000 feet that he was 
2.2 degrees shallow of his desired flight path. In the 
second case h.e discovered that he 1.9 degrees steep of the 
desired track. For each dive angle error two alternative 
techniques were used to re-establish the desired 40° track 
over the ground. The first was to shallow or steepen 
further an amount which would allow re-establishment of the 
desired track immediately prior to the release altitude of 
5,000 feet. The second technique was to shallow or steepen 
more than the first to allow re-establishment on the desired 
flight path early enough for the major protion of the air- 
frame short period oscillations to die out prior to the 
release altitude of 5,000 feet. Figures 4 through 7 illus- 
trate the two error cases and the flight paths used to cor- 
rect for these errors. 

B. METHOD OF AIRCRAFT CONTROL 

In shallowing or steepening the aircraft dive angle, 
the full system of nine equations was used to solve for the 
time response of u, a, 0, q and y to the commanded perturba- 
tion dive angle y^. When the correction dive angle or the 
re-established 40° dive angle was approached, the reduced 
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Figure 4. Run 1, 2.2° Shallow, 5° Correction. 
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Horizontal Range from Release Point (feet) 
Figure 7. Run 4, 1.9° Steep, 8° Correction. 



system consisting of only the four airframe equations was 
used to solve for the perturbation quantities. This was 
done for two reasons: 

1. The simplified transfer functions of the longitudi- 
nal control system were invalid in the pitch rate range be- 
tween -0.5® per second and 0.5® per second. Obviously a 
constant dive angle or zero pitch rate falls within this 
range. 

2. Since it was the short period angle of attack per- 
turbation which was of interest, the radian frequency for 
this mode of oscillation was such that airframe damping 
rather than pilot control movement is the factor which 
dampens the oscillation. ' 

The pilot then flew the aircraft from the initial dive 
angle to the correction dive angle through the control sys- 
tem, anticipated achieving the new desired dive angle based 
on pitch rate and set the stabilator to the reference posi- 
tion (perturbation 6g = 0) for the new dive angle. 

After the aircraft had flown the required time on the 
adjusted flight path to approach the correct 40® dive angle, 
the controlling sequence was repeated with the pilot re- 
establishing the 40® dive and the airframe dynamics in the 
new reference state maintaining it. 

C. COMPUTER THEORY FOR AIRCRAFT CONTROL 

Reference 6 contains computer programs for computational 
assistance in the study of linear control theory. The 
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programs used were GTRESP, CALCU, RUNGE, TRESP, YDOT, Y8VSX 
and user supplied subroutine RFIND. Program listings and 
explanation are contained in Appendix D. The purpose of the 
programs is to furnish the graphical time response to the 
following system: 

{x(t)} = [A] {x(t)> + {b} [u(t)] 

[u(t)] = K(r(t) - {k}”^ {x(t)}) 
y(t) = {c}^ {x(t)}. 

In the case under consideration, {x} is the vector 
of perturbation quantities, [A] is the nine-by-nine or 
four-by-four plant matrix defining the airframe and control 
system, {b} is the control vector, K = is the pilot stick 
force per radian of difference between actual and commanded 
di ve angle, r is the commanded perturbation dive angle, {k} 
is the vector of feedback coefficients, and y is the actual 
perturbation dive angle made up of the output vector {c} 
multiplied by the airframe perturbation vector 

(Y = {0 -1 1 OHu a 0 

Following is the detailed breakdown of one of the 
situations considered. At time zero the pilot at 9,000 
feet altitude in a 40° dive determined that he was 885 feet 
horizontal range short of the correct position for that al- 
titude in the delivery maneuver. This was equivalent to 
being 2.2° shallow. He elected to correct for this by shal- 

lowing 5°, flying at a 35° dive angle until approaching the 
> 

desired track, steepening back to a 40° dive angle and flying 

at that 40° angle to release. 
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GTRESP was called with the following data: 

1. the nine-by-nine plant matrix with airframe data 
referenced to a 40® dive angle and the control system, 

2. the control vector, 

3. the output vector, 

4. the feedback coefficient vector, 

5. Kp based on required dive angle change of 5°, 

6. equal to +5° in radian measure, and 

7. zero initial perturbation vector referenced to 
40° dive conditions. 

GTRESP then used the other subroutines to perform a 
fourth order Runge-Kutta integration to compute and plot the 
time response of the nine perturbation quantities. 

At 0.975 seconds the perturbation dive angle was 
approaching 5° or T was approaching 35°. The airframe per- 
turbation quantities for time equal to 0.975 seconds were 
transformed to perturbation quantities for a 35° dive ref- 
erence state, and GTRESP was called with the following in- 
put data: 

1. the reduced plant matrix of airframe data only 
referenced to a 35° dive, 

2. the output vector, and 

3. aircraft perturbation quantities referenced to 
35° dive at time equal to 0.975 as the initial condition 
vector. 

Failure to input the control vector, feedback coeffi- 
cients, K and bypasses the control system which is not 
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accurate for zero and very small pitch rates and simulates 
the pilot controlling the aircraft by setting and maintain- 
ing the zero perturbation stabilator position for the 35* 
dive reference condition. Again GTRESP computes and plots 
the time response of the aircraft perturbation quantities. 

At 7.2 seconds total elapsed time the aircraft was 
approaching the desired '40° track over the ground and GTRESP 
was called with the large plant matrix for control back to 
a 40° dive. Input was the same as for the first maneuver 
except that the aircraft dynamics portion of the nine-by- 
nine matrix was based on the 35° dive since that was the 
condition from which perturbations occurred, and the ini- 
tial condition perturbations remained referenced to the 35° 
dive . 

GTRESP computed and plotted the perturbation quan- 
tities through the transition back to a 40° dive. At 8.175 
seconds total elapsed time the 40° dive on the correct 
flight path was re-established and it was simulated that the 
pilot neutralized longitudinal control by calling GTRESP, 
with the reduced plant matrix and initial perturbations 
transformed to those for a 40° dive reference state, to 
calculate and plot the perturbation dive angle, velocity, 
angle of attack, pitch angle, and pitch rate to release 
point . 

On this particular run it was arranged to have the 
aircraft re-achieve the desired 40° dive angle at exactly 
the 5,000 feet release altitude. Therefore, the angle of 
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attack perturbation, converted to mils, existing at 8.175 
seconds is in fact the release sight picture error that the 
pilot would have had on this maneuver at release altitude. 
Pertinent aircraft position and gunsight error information 
concerning each of the four bomb delivery maneuvers simu- 
lated is listed in Table II. 

Computer output showing time history in 'tabular form 
of aircraft perturbation during the four phases of each of 
the four runs tested and in graphical form for the final 
delivery phase is enclosed in the Computer Output Section. 

D. RESULTS OF THE FOUR RUNS 

Table II shows that the total time available for ma- 
neuvering after detecting a dive angle error, 4,000 feet 
above release altitude, is very small for a 40° 500-knot 
dive bombing run. Therefore, any corrections made must be 
both timely and correct, if they are to succeed in increas- 
ing the accuracy of the bomb impact. Table II also shows 
that by merely increasing the correction dive angle from 
5° to 8°, and thereby re-establishing the desired dive path 
1.75 seconds prior to release rather than at release for 
the shallow run and 1.4 seconds prior to release rather than 
at release in the steep run, reduces the gunsight position 
error from 83 to 9 mils for the shallow case and from 80 to 
3 mils for the steep case. 

Examination of the graphical time histories of perturba- 
tion angle of attack, presented in Appendix C, shows that 
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Sight 

Height Horizontal Horizontal Total Error 

RUN 1 Alt. Above Distance Distance Elapsed at 

2.2° Shallo\v Lost Release Travelled to Release Time Release 

5° Correction (ft) (ft) (ft) (ft) (sec) (mils) 
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Table II. Aircraft Position During Dive Angle Corrections and Gunsight Error at Release Alt. 



any additional time spent on the desired dive path rapidly 
decreases the angle of attack perturbation, thereby rapidly 
increasing the accuracy of the gunsight position relative 
to actual bomb impact point. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



Table II shows the elapsed time in each of the four 
weapons delivery maneuvers tested at which the aircraft 
finally becomes established in a 40® dive on the proper 
track over the ground for an accurate bomb hit. The Com- 
puter Output shows the FRL perturbation angle of attack at 
numerous time increments in each run. It is obvious for 
both the initially steep and initially shallow runs that 
the 5° correction dive angle, although sufficient to correct 
the aircraft back to desired dive path by release altitude, 
is not large enough to correct back to the desired dive 
path soon enough to allow for sight picture settling time, 
due to angle of attack oscillations of the short period 
longitudinal mode. However, by merely Increasing slightly 
the magnitude of the correction dive angle change, it was 
possible to reposition the aircraft on the desired flight 
path early enough before the release point to allow time for 
the sight picture to stabilize sufficiently to achieve an 
acceptable bomb impact. Examination of the graphs showing 
time history of angle of attack perturbations during the 
final phase of the delivery maneuver shows that even earlier 
establishment on desired dive path is desirable and would 
increase the accuracy of the sight picture further. 

In actual practice, a pilot on detecting a dive angle 
error would not consciously attempt to hold a dive angle 
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correction for a specific amount of time with the idea that 
in so doing he would re-establish his aircraft on the proper 
flight path. If he could do this, there would be no need 
for the gunsight. All the cases described herein would 
actually achieve very accurate bomb impacts, since for an 
aircraft at the correct velocity and in the desired dive 
angle along the proper path over the ground, very large 
angle of attack perturbations produce bomb impact errors 
only insofar as they change the aircraft dive angle. The 
error which would actually be produced in the case where 
all delivery parameters are met except for an angle of at- 
tack error would be that the gunsight reticle would be far 
off the target. 

However, since the pilot in a visual dive bomb delivery 
uses visually acquired data for dive angle and aircraft 
positioning, he is entirely dependent on his gunsight for 
accurate aircraft positioning information; and during a 
period when he is receiving inaccurate gunsight information, 
he has no good indication of his position relative to the 
desired release position. 

It has long been standard practice in visual air-to-air 
gunnery with a gyro-stabilized gunsight, which develops its 
own lead angle based on the normal acceleration of the fir- 
ing aircraft required to track the target, to teach students 
not to squeeze the trigger until they have maintained a 
smooth tracking solution for two seconds on the target. This 
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was to allow angle o£ attack oscillations from a changing 
pitch rate to die out prior to firing. 

The need to apply this same technique to visual air-to- 
ground delivery is just as important, or even more so when 
it is considered that a gyro-stabilized sight has inertia 
to decrease the magnitude of the effect' of aircraft angle 
of attack perturbations; whereas the sight line of a mech- 
anically fixed sight moves directly with angle of attack 
perturbations an amount equal to the movement of the FRL. 

Nor does a pilot in air-to-ground delivery have the' option 
to wait until perturbations die out if this will take him 
below the planned release altitude, since the computed sight 
angle is good for only the planned delivery parameters. 

What he must do is to detect dive angle errors sufficiently 
early in the run and make corrections sufficiently large 
that all corrections are completed and all significant 
perturbations eliminated prior to the planned release height. 

This study pertained specifically to the F-4 aircraft in 
a 40° dive, MK 82 LDGP Bomb, delivery. However, since all 
Navy fighter and light attack jet aircraft have similar 
short period characteristics with slight variance in damping 
ratio and oscillation frequency, and since all fixed-sight 
air-to-ground weapons delivery maneuvers have similar error 
sensitivities (differing only in the magnitude of the range 
error for each mil of aimpoint error) , the results are ap- 
plicable to all fixed gunsight, visual, air-to-ground weapons 
deliveries in all Navy tactical aircraft. As such, it is 
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recommended that similar studies be made to encompass all 
air-to-ground weapons deliveries from Navy tactical aircraft 
and the results be used in the lecture portions of air-to- 
ground training in both fleet and replacement squadrons. In 
addition a graphical display of the results herein can be 
easily demonstrated to a pilot in flight in the following 
manner. Trim the aircraft for steady level flight with the 
gunsight set on the distant horizon. Make a small but sharp 
longitudinal displacement of the control stick and return it 
to neutral. Notice that, although the change in aircraft 
track (speed and altitude) is extremely small, the magnitude 
of the gunsight oscillation above and below the horizon is 
significant; and that while the time to damp is short, it is 
long enough to lose several hundred feet, if it were in a 
40° 500-knot dive. 

Awareness of this phenomenon can do much to improve the 
proficiency of the not atypical pilot who, on being ques- 
tioned about pipper placement after an unsuccessful air-to- 
ground weapons training flight, is not really able to tell 
an instructor whc.o his pipper was positioned at release on 
most of his runs, this because of late corrections resulting 
in receipt of meaningless data from the gunsight. 
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APPENDIX A: DERIVATION OF LINEAR CONTROL FEEDBACK COMPUTER 

REPRESENTATION OF F-4 LONGITUDINAL CONTROL 
SYSTEM 



The longitudinal control system of the F-4 aircraft is 
shown below in diagrammatical form, 




where , 



Feel System = G2(S)=27.1[ 



S"+ 5.637S+565. 5 



1 f inches stick travel s. 
^ pound force ^ 



Control Stick Gearing s G,(S)=-0. 05864 ( radians stabllator travel 

^ inch stick travel ^ 



Stabilator Actuator = G.(S)= 20 [ g— ^ ^ ] 

Bobweight = H,(S)= 0.1663 ( P°unds bobweight force j 



radians stabilator travel. 



stab. Aug. System . H^(S)= 0.15 [ 1 ( ^d^Wser pUch r^^r > 



With the addition of airframe dynamics and a pilot, modeled 



by Kp pounds of stick force per radian of error becween actual 



and desired pitch angle, the representation becomes 
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where the additional transfer functions are. 



Pilot force/radian e. = G. (S) = K 

b 1 P 



( pounds force ) 
radian 



g 

a^/^s ti'ansfer function 



= H^(S) = 861. 5[ 4 



S^+1.213S^+27.98S^+2.272S+0.04339 



S '" + 2. 8025^+15. 19S^ + 1. 179S + 0. 0 224 



]( 



g 



radian 



) 



0/6g transfer function 



Gr(S) = 



- 24 . 64[-4 



S^+1. 2S+5 , 367 



S'^ + 2.802S^+15. 19S^ + 1.179S+0.0224 



.1 (- I’adian . 
■* ^radian 



q/0 transfer function = H^CS) = S 



( 



radian/ sec 
radian 



) 



• A block diagram algebra reduction of the control system 
as described in Reference 5 follows: 




44 







f 7 » 




I 



E2 = El - HiH^Sg 

' ®2‘^2®3'’4 “ ■ ‘'2''3‘^4^1^2*S 



^ ■ °2°3°4 

1 1 * GiGjG^HiHi 




E 



1 




'=sH3»4 , 

G 2 G 3 «S 



^S = 



SS<=4 

1 . GiGiG^HiHi 



^2^3^4 

«s ,. 

Fjt - ^ GiGjGiGjHjH^ 

• ■ GjGsU-GiGjG^HiHi) 

^S ^ ^2^3^4 

Fst ^ ■" ^2^3^4^1^2 ‘ ^4^5^3^4 
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e 



^ST ■ ®e^i 

f _ c 

~ ^ST 1 + G.G_G.H,H. - G.G.H.H. 

23412 4534 

S ' 1 ♦ G^GjG^HjH^- G^G^HjH^ 

Cl * GjGjG^H^Hj - G^GjHjH^) 63 ( 8 ) = G^G^GjG^ [ ( 6 ^- 6 ) ( 8 ) ] 



In this final form the control system in the Laplace 
domain can be readily transformed to its state space repre- 
sentation of a linear differential equation for- 6 g as a 
function of time and (9^ - 0 ) . However reducing all the 
transfer functions which are polynomial fractions in S to a 
least common denominator would result in a forty second 
order linear differential equation. Therefore, representa- 
tion of some of the transfer functions that would yield a 
resulting differential equation more amenable to numerical 
solution but still accurate, was sought. 

0 

— (S) s G^(S) was not disturbed since it represented 
S 

actual airframe dynamics which was under investigation. 
Neither was the Stab. Aug. System = H^(S). However, a fre- 
quency response investigation of the control stick h 02 ( 8 ), 

g 

the stabilator actuator = G^(S) and a^ = Hj^(S), using the 
FRESP computer program of Reference 6 , yielded Bode Diagrams 
of transfer function magnitudes and phase angles showing 
that over the pitch rate range of interest ( 0.01 to 0 . 2 ) 
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radians/second - a normal acceleration range up to ± 5 . 2 ) the 
above transfer functions could accurately be represented by 
constant magnitudes with zero phase angle. 

With these simplifications, the transfer functions are 



as follows : 






G^( 8 ) = 


K 

P 


H^( 8 ) = 1600 


G^CS) = 


0.04792 


H 2 ( 8 ) = 0.1663 


GjCS) = 


0.05864 


H 3 ( 8 ) = 8 


G4(S) = 


1.0 


H 4 ( 8 ) = 0.15[ ^ 



G^(S) = -24.64 [ 



+ 1.25 + 5.367 

+ 2.8025-^ + 15.195^ + 1.179S + 0.0224 



The Laplace domain representation for 6 g as a function 
of time and (6 - 0 ) becomes 

(S^ + 18.455^ + 35.578^ + 95.08^ + 1.2018 + 0.0224) 65 CS) 

- 0.01113K (S^ + 3.8028^ + 17.998^ + 16.378^ + 1.2018 
P 

+ 0.0224) [0^ - 0 ( 8 )] . 

The value of Kp, pilot force per radian of pitch angle 
error, was varied depending on what gain was required to 
achieve an average pitch rate of 0.115 radians per second 
which is equivalent to a change in normal acceleration from 
steady state conditions of 3g (typical in fact of the pitch 
rate a pilot would use in making a dive angle correction) . 
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In addition, the pilot was assumed to be aware that 
actual flight path angle (dive angle) lagged the pitch at- 
titude angle during a changing pitch rate by the amount of 
instantaneous angle of attack. Since the purpose was dive 
angle change not just pitch attitude change, 0-a or y was 
fed back in the control system to solve for 6^ as a function 
of time and (y^-y) rather than (6^-0). This is the repre- 
sentation shown in the matrix feedback differential equation 
of Figure 3. 

The simplification of the transfer functions listed 
above necessitated using airframe dynamics alone, in a 
simulated stick fixed condition, rather than coupled with 
the control system after the desired dive angle was approach- 
ed, since very small pitch rates were outside the limits for 
which the simplified forms of (S) , G^(S) and H^(S) applied. 
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APPENDIX B: NUMERICAL VALUES OF CONTROL VECTOR AND ELEMENTS 

OF THE PLANT MATRIX FOR REFERENCE DIVE ANGLES 
USED 

The elements of the control vector derived in Sections 
II-B and II-C are invariant with dive angle. Values of 
elements of the control vector are listed in Figure 8. 

Derivation of formulae for the elements of the plant 
matrix is in Sections II-B and II-C, and the matrix is 
shown in symbol form in Figure 3. Definitions of the ele- 
ments is in "Definitions of Symbols" at the front of this 
study. Calculation of the airframe dynamics entries was 
accomplished in the following manner: 

1. Dimensionless stability axis data were obtained 
from References 1 and 2. 

2. Dimensionless body axis figures were computed from 
the "Definitions of Symbols", using stability axis data 
and reference state angle of attack. 

3. Dimensional body axis figures were calculated using 
dimensionless body axis figures, dimensional derivative 
definitions, and applicable reference state data. Tables 
III and IV are compilations of the unperturbed 40° reference 
state conditions, and aerodynamic coefficients/stability 
derivatives for each reference dive angle used. Figure 8 
lists the values of the plant matrix elements for each dive 
angle . 
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Since the aircraft was assumed in a constant speed 
(except for perturbations) 500 knot dive, rather than ac- 
celerating from 430 KTAS to 500 KTAS, lift and drag coef- 
ficients from References 1 and 2 were invalid. These 
coefficients were calculated for each dive angle from free 
body force balance diagrams. 



» 
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Table III. 40° Dive Reference Conditions and Stability Axis 

Aerodynamics Coefficients and Stability Derivatives 

To , = -40° = -0.6981 radians 

Vj = 844.4 ft/sec 

0 

Uq = 844.4 ft/sec 

Wq = 3.369 ft/sec 

Mo = .775 

©0 = -39.77° = -0.6941 radians 

a = 3.99 mils = 0.2286° = 0.00399 radians 

0 

m = 1210 slugs 

W = ^ 38920 lb 

p = 0.001933 slugs/ft^ 

S = 530 ft^ 

c = 16.04 ft 

q = 689.1 Ib/ft^ 

g = 32.17 ft/sec^ 

I = 122200 slug ft^ 

y 

c, = 5.25° 

0 

1^, = -0.57 ft 

th 

Ig = 39.4 ft 

Cg = 0.06535 

Sa = 0.1464 

^DM ^ 0 

» 

Cg = 0.1391 
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3.36 



^La 

^La 

^LM 

^Lq 

<=L63 

Cm 

c 

ma 

C • 

ma 

^mM 

C 

mq 

Cm6s 



0.7 

0 

2.22 

0.355 

0 

-0.287 

- 1.0 

-0.0314 

-2.37 

-0.515 
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Table IV. Body Axis Force Coefficients and Stability Deriva- 
tives . 



r 


-32° 


-35° 


-40° 


-45° 


0 

00 

1 


Cx 


0.05026 


0.0611 


0.06479 


- 0.07144 


0.07717 


^Xa 


-0.002913 


-0.002064 


-0.006367 


-0.006422 


-0.01468 


^XM 


0 


0 


0 


0 


0 




-0.001416 


-0.001416 


-0.001416 


-0.001416 


-0.001416 


s 


0.1259 


0.1317 


0.1394 


0.1495 


0.1628 


Sa 


3.411 


3.422 


3.425 


3.432 


3.437 


Sa 


0.7 


0.7 


0.7 


0.7 


0.7 


^NM 


0 


0 


0 


0 


0 


^Nq 


2.22 


2 .22 


2 .22 


2.22 


2.22 


S6s 


0.355 


0.355 


0.355 


0.355 


0.355 


X 

u 


-0.03541 


-0.04345 


-0.04633 


-0.05129 


-0.05566 


X 

a 


0.6541 


0.4344 


1.766 


1.911 


4.534 


X 

q 


-0.0254 


-0.0254 


-0.0254 


-0.0254 


-0.0254 


X6s 


0.4274 


0.4274 


0.4274 


0.4274 


0.4274 


z 

u 


-0.08029 


-0.08759 


-0.09478 


-0.1061 


-0.118 


z 

a 


- 1010. 


- 1022. 


- 1034. 


- 1045. 


- 1056. 


Z- 

a 


-2.007 


-2.007 


-2.007 


-2.007 


-2.007 


Z 

q 


-6.365 


-6.365 


-6.365 


-6.365 


-6.365 




-107.2 


-107.2 


-107.2 


-107.2 


-107.2 


M 

u 


-0.001317 


-0.001317 


-0.001317 


-0.001317 


-0.001317 


M 

a 


-13.76 


-13.76 


-13.76 


-13.76 


-13.76 


M* 

a 


-0.4553 


-0.4553 


-0.4553 


-0.4553 


-0.4553 


M 


-1.079 


-1.079 


-1.079 


-1.079 


-1.079 



q 
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M6s 


-24.69 


-24.69 


-24.69 


9T 

3M 


0 


0 


0 


T 

u 


0 


0 


0 



-24.69 -24.69 

0 0 

0 0 
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PLANT MATRICES FOR EACH REFERENCE DIVE ANGLE AND THE CONTROL VECTOR 
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FIGURE 8 (CONT 
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TIME RESPONSE OF PERTURBATION QUANTITIES DURING PUSH OVER FROM 35 TC 40 DEGREES 
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TIKE RESPONSE OF PERTURBATION QUANTITIES CURING 40 DEGREE DIVE 
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GRAPHICAL TIME RESPONSE OF PERTURBATION GAMMA AND ALPHA AFTER RE-ESTABLISHMENT IN DESIRED 40® DIVE 
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GRAPHICAL TIME RESPONSE OF PERTURBATION GAMMA AND ALPHA AFTER RE-ESTABLISHMENT IN DESIRED 40® DIVE 

VARIABLE SYMBOL 
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PUN 3, 1.9 DEGREES STEEP, 5 DEGREE CORRECTION 
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TIME RESPONSE OF PERTURBATION QUANTITIES CURING 45 DEGREE DIVE 
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TIME RESPONSE OF PERTURBATION QUANTITIES DURING PULL UP FROM 45 TO 40 DEGREES 
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TIME RESPONSE OF PERTURBATION QUANTITIES CURING 40 DEGREE DIVE 
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GRAPHICAL TIME RESPONSE OF PERTURBATION GAMMA AND ALPHA AFTER RE-ESTABLISHMENT IN DESIRED 40“ DIVE 
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TIME RESPONSE OF PERTURBATION QUANTITIES CURING 48 DEGREE DIVE 
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TIME RESPONSE OF PERTURBATION QUANTITIES DURING PULL UP FROM 48 TO 40 DEGREES 
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TIME RESPONSE OF PERTURBATION QUANTITIES CURING 40 DEGREE DIVE 
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COMPUTER PROGRAMS 



USER MAIK PROGRAM, EMPLOYED TO START 
EACH PHASE OF A DELIVERY MANEUVER 



CALL GTRESP 
CALL GTRESP 
CALL GTRESP 
CALL GTRESP 
01 STOP 
END 



TYPICAL USER SUBROUTINE, EMPLOYED TO SET 
START TIME FOR EACH PHASE OF MANEUVER 



SUBROUTINE RFIND(T,R) 
IFCT.GT .0.975J GO TO 2 
R=. 06725 
RETURN 

2 IF(T.GT .7.2) GO TO 3 
R = 0. 

RETURN 

3 I F(T.GT .8 . 175) GO TO 4 
R=-0. 08725 

RETURN 

4 IFIT.GT .15.675 ) GO TO 5 

R=0. 

5 RETURN 
END 
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PRIMARY SUBROUTINE, EMPLOYED TO CALL ALL OTHER ROUTINES 
FOR CCMPUTING, BY FOURTH ORDER RUNGE KUTTA NUMERICAL 
INTEGRATION, AND PLOTTING THE TIME RESPONSE CF THE 
PERTURBATION QUANTITIES TO A COMMANDED DIVE ANGLE 

CHANGE 



SUBROUTINE GTRESP 
GRAPHICAL TIME RESPONSE (GTRESP) 

SUBPROGRAMS USED- CALCU, RUNGE, TRESP, YDOT 
INTEGER CHARdS) 

CCMMON I PLOT, I VAR ( 10) 

DIMENSION A( 10, 10) ,C( 10) , B( 10) , AK ( 10 ) , X (1 0 ) , NAME ( 5 ) 
DATA CHAR ( 1) , CHAR ( 2 ) , CHAR ( 3 ) , CH AR (4) ,CHAR( 5 ) , 

=^CHAR(6) ,CHAR(7) ,CHAR(8) ,CHAR(9) ,CHAR( 10), CHAR (11), 
«CHAR(12),CHAR(13) ,CHAR(14),CHAR(15)/2H 1 , 2H 2,2H 3,2H0 
*4,2H 5,2H 6,2H 7,2H 8,2H 9,2H10,2H E,2H U,2H Y,2H R, 
*2H / 

3 FORMAT (6F10.3) 

1000 FORMAT (IHO, lOX, 8HTZER0 = ♦ F10.6,10X, 5HTF = ,F10.6 
*/ 

* llX,5hDT = ,F10.6,13X,7HFREQ = 15) 

1001 FCRMaT (1H0,10X,13H the a MATRIX /) 

1002 FORMAT (6( 1PE20.8) ) 

1003 FORMAT (1H0,10X,1VH INITIAL CONDITIONS /) 

1004 FCRMAT (1H0,10X,13H THE B MATRIX /) 

1005 FORMAT (1HO,10X,16H FEEDBACK COEFF. /) 

1006 FORMAT (iH0,10X,8H GAIN = , 1PE20.8 ) 

1007 FORMAT (1H0,10X,13H THE C MATRIX /) 

1008 FCRMAT (8A2) 

1009 FCRMAT( 5X , 25HPR0BLEM IDENTIFICATION - ,5A4) 

1010 FORMAT! IHl ,4X,23HGRAPHICAL TIME RESPONSE) 

1011 FORMAT! /5X,45( IH'T) ) 

10READ(5,1) (NAME! I ), 1 = 1,5) ,N 

1 FORMAT ( 5A4, 12) 

DO 60 1=1,8 

60 I VAR( I J =CHAR(15) 

PRINT 1010 

PRINT 1009, (NAME( I ), 1=1,5) 

PRINT lOil 

PRINT 1001 

DO 2 I = 1, N 

READ 3, (A( I , J) , J=1,N) 

PRINT 1002, (A(I,J), J=1,N) 

2 CONTINUE 

READ 3, ( B( I ) , 1=1 ,N) 

PRINT 1004 

PRINT 1002, (B(I),I=1,N) 
read 3, (C( I ) , 1=1 ,N) 

PRINT 1007 

PRINT 1002, (C( I) , I=1,N) 

READ 3, ( AK( I ) , 1=1, N) 

PRINT 1005 

PRINT 1002, (AK( I ) , 1=1 ,N) 

READ 3, GAIN 
PRINT 10C6, GAIN 
READ 3, (X(I) ,I=1,N) 

PRINT 1003 

PRINT 1002, (X(I),I=1,N) 

READ 3, TZERG, TF,DT,FREQ 
IFG=FREQ 

PRINT 1000, TZERG, TF ,DT,IFQ 
PRINT 1011 

READ lOOS, ( IVAR( I ) ,I =1,8) 

DC 40 1=1,8 
DO 30 J=1 , 15 

IF (IVAR(I)-CHAF(J)) 30,25,30 
25 iVAR(I)=J 
GO TG 40 
30 CONTINUE 
40 CONTINUE 
MIN = 1 
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MAX=8 
M = 8 

AIS DO 42 I=MIN,MAX 

IFdVAKiD.NE.iS) GO TO 42 
M=MAX-1 

IFU.GT.M) GO TO 42 
DC 43 J = I f M 

43 IVAR( J)=IVAR( J+1) 

GO TO 431 

42 CONTINUE 
GO TO 432 

431 MIN=I 
MAX=M 

GO TO 419 

432 I PLOT = M 
IF(IPLOT.LT.2) GO TO 50 
LIM=1 PLOT-1 

DO -44 -L=l, LIM 
MIN=I+1 

DO 44 J=MIN, I PLOT 
I F( IVARd )-IVAR( J ) ) 44,44,45 
45 IHGLD=IVARd) 

IVAR{ I ) =IVAR( J ) 

I VAR( J) =IHOLD 

44 CONTINUE 

50 CALL TFESPI A,X,B,AK,TZERO,TF,OT, IFQ,N,GAIN,C) 
20 CONTINUE 
RETURN 
END 
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SUBROUTINE TRE SP ( A t Y , B , AK , X , XMAX , H , I FREQ , N , GA IN , C ) 
THIS SUBROUTINE COMPUTES AND PLOTS TIME RESPONSE 
USING CALCU, PURGE, YOCT AND Y8VSX 
INTEGER CHARI 15) 

COMMON IPLOT, IVAR( 10) 

DIMENSION SKJ (101,9) ,C(10) 

DIMENSION FN(iO), Y( 10) , A ( 10, 10) , B ( 10 ) , AK( 10 ) 

DATA CHAR( 1) , CHAP ( 2 ) ,CHAR (3 ) , CHAR ( 4) , CHAR ( 5 ) , 
^CHAR(6) ,CHAR(7 ) , CHAR ( 8 ) , CHAR(9 ) , CHAR ( 10) ,CHAR( 11) , 
=C=CHAk(12),ChARll3),CHAR( 14) , CHAR ( 15 ) / IHl , 1H2,1H3,1H4, 
’('IhStlHb , 1H7, IHb, 1H9 ,1HA, IHE, IHU, IHY, IHR , IH / 

24 FORMAT ( 2F 10 . 0 , 2 1 1 0 ) 

25 FORMAT (8F10.0) 

28 FORMAT ( / / , 8X , IHT , 12X, 4HY ( T ) , lOX , 4HU( T ) , 4X , 

* 6(5X, IHX, I1,4H(T) ,3X)) 

29 FORMAT ( 9( 1PE14.6 ) ) 

1000 FORMAT! /, 5X,33HMAX1MUM NUMBER OF POINTS EXCEEDED /» 
PRINT 28, ( J,J = 1,N) 

I 1=0 
J=0 

KCUNT=I FREQ 

300 CALL CALCU(Y,U,X,N, AK,GAIN,R,0) 

KCUKT = KOUNT + 1 
IF (KOUNT - IFREQ) 50, 350, 350 
350 KOUNT = 0 

450 P1=0.0 

DO 451 1=1, N 

451 P1=P1+C ( I )*Y( I ) 

PRINT 29, X,P1,U, ( Y(M) ,M=1,(N) 

WRITE (7,31)X,P1, (Y(M) ,M=1,4) 

31 F0RMAT(6( 1PE12.3) ) 

IF (IPLOT. EQ.O) GO TO 21 
J=J+1 

IF(J.GT.lOl) GO TO 222 
SKJ( J,1 )=X 
DO 40 1 = 1 , IPLOT 
MM=IVAR(I J 

IF(MM.EQ.O) GO TO 40 
IF(MM.GT.iO) GO TO 35 
SKJ( J, I+l ) =Y(MM) 

GO TO 40 

35 KN0W=MM-10 

GO TO (36,37,33,39), KNOW 

36 SKJ(J, I+1)=R-P1 
GO TO 40 

37 SKJ(J, !+i) =U 
GO TO 40 

38 SKJ(J, 1+1)=P1 
GO TO 40 

39 SKJ(J, I+1)=R 

40 CONTINUE 
21 CONTINUE 

50 CALL RUNGE (N, FN, H, X, Y, L,II) 

IF(L-l) 100, 200, 100 

200 CALL C ALCU(Y,U,X,N,AK,GAIN,R,1) 

CALL YDOT(A,Y,FN,B,U,N) 

550 GO TO 50 
222 PRINT 1000 
GO TO 400 

100 IF(X-XMAX) 300,300,400 
400 I F( IPLOT. EQ.O ) GO TO 403 
PRINT oOO 

600 FORMAT ( 1 H 1 , 50X , 1 5HS YSTEM RESPONSE//) 

PRINT 601 

601 FORMAT ( 48X , bH VAR I ABL E , 8X, 6HSYM BOL// ) 

DO 606 1 = 1, IPLOT 

MM = IVAR(I ) 

IF(MM.GT.IO) GO TO 603 
PRINT 602, I VARII ) ,CHAP.(MM) 

602 FORMAT ( 5 1 X , IHX , I 2 , 13X , A1 ) 

60 TO 608 

603 KK=1VAR(I )-lO 
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GG 70(604,605,606,607), KK 

604 PRINT 610 
GO TO 608 

605 PRINT 611 
GG TO 6 08 

606 PRINT 612 
GC TO 603 

607 PRINT 613 

608 CONTINUE 

610 FORMAT ( SOX , 5HERR0R , 12X , IHE ) 

611 FORMAT (49X, 7HC0NTR0L, ilX, IHU) 

612 FORMAT ( 50X, 6HGUTPUT , IIX, IHY) 

613 FORMAT { SOX , 5H i N PUT , 12X , IHR ) 
CALL Y8VSX(SKJ, J, IPLCT,10) 

403 RETURN 
END 
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SUBROUTINE C ALCU ( X t U , T , N, AK t GA I N , R , J ) 

THIS SUBROUTINE COMPUTES THE REFERENCE AND CONTROL 
INPUTS 

DIMENSION X(10) , AK(IO) 



CALL RFINO(T,R) 



END OF ROUTINE TO DEFINE R(T) 
U = R 

00 1 I = It N 
1 U = U - AK( n«X( I ) 

U = U«GAIN 

RETURN 

END 



SUBROUTINE RUNGE (N,FN, H, X» Y, L,I) 

C FOURTH ORDER RUNGE KUTTA INTEGRATION ROUTINE 

DIMENSION Y(600)j SAVEY(600), PHI(600), FN(8) 
1 = 1 + 1 

GO TO ( 1, 2, 3, 4, 5) , I 

1 L = 1 
RETURN 

2 DO 600 J=1,N 
SAVEY(J) = Y(J) 

PHKJ) = FN(J) 

600 Y(J) = SAVEY(J) + .5-H^FNU) 

X = X + .5*H 
L = 1 

RETURN 

3 DO 700 J=1,N 

PHI(J) = PHKJ) + 2.*FN(J) 

700 Y(J) = SAVEY(J) + . 5=i'H=4=FN ( J ) 

L = 1 
RETURN 

4 DC 800 J=1,N 

PHKJ) = PHKJ) + 2.*FN(J) 

800 YU) = SAVEY(J) + H*FN(J) 

X = X + .5*H 
L = 1 
RETURN 

5 00 900 J=ltN 

900 Y(J) = SAVEY(J) + ( H/6 , ) «( PH I ( J ) + FN ( J ) ) 

L = 2 
I = 0 

return 

END 
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SUeROUT IME YDOT( A,Y,XDOT,BtU,N) 

C THIS SI.IBROUTIME 13 USED TO COMPUTE DERIVATIVES FOR 

C RUNGE 

DIMENSION Y(10), AU0,1G), B(10), XDOT(IO) 

DC 2 I = 1, N 
XCGT( I) = 0. 

DO 1 J = 1, M 

XDCT(I) = XDOT(I) + A(I,J)«Y(J) 

1 CONTINUE 

XDOT( I ) = XOOT( I ) + B( I )*U 

2 CONTINUE 
RETURN 
END 



SUBROUTINE Y8VSX ( A , N, M, NGRI D) 

C THIS SUBROUTINE PLOTS UP TO 3 VARIABLES VERSUS TIME 
C THIS IS BOTH A X-T AND X-Y PLOT ROUTINE 
COMMCN IPLOT, I VAR( 10) 

INTEGER CHAR(15) 

DIMENSION A( 101 ,9) , ABSCA( ll)tKAXIS(lOl), CRDINdl ) , 
«TEMPYI9 ) 

DATA CHAR m , CHAR (2) ,CHAR ( 3 ) , CHAR (4) , CHARI 5 ) , 
«CHAR(6) , CHAR (7) ,CHAR(8) , CHARI 9) , CHAR I 10) ,CHARI11) , 
*CHARI 12) , CHAR I 13 ) , CHAR I 14) , CHARI IB ) / 1 HI , 1 H2 , 1H3 , 1H4 , 
*1H5,1H6, 1H7, 1H6, iH9tlHA,lHE ,1HU,1HY,1HR,1H / 

DATA ISTAR, 1 I? IPER, IDASH, IBLANK/lh=!=, IHI , IH. , 1H-, IH / 
ICO FORMAT I////,9X, 11I1PE10.2)) 

101 FORMATI 1PE13.2,2X, lOlAl) 

102 FCRHATI 15X, iOlAl ) 

YMAX=AI 1,2) 

YMIN=AI 1,2) 

MP1=M+1 

DO 4 J=2,MP1 
DO 4 1=1, N 

IFIYMAX-AI I, J) ) 1,2,2 

1 YHAX=AI I, J) 

2 I FIYMIN-AI I, J) ) 4,4,3 

3 YMIN=AI I, J) 

4 CONTINUE 

5 YSHFT = YMIN=!=100.0/ I YMAX-YMIN) 

6 NM1=N-1 

DO 8 1=1, NMl 

iPl=I+l 

Du 8 K=IP1,N 

IFIAIK, D-AI 1,1)) 7,8,8 

7 DC 88 J=1,MP1 
ATEMP=AI I , J) 

AI I ,J) =AI K, J) 

AIK, J)=ATEMP 

88 CONTINUE 

8 CONTINUE 
XMIN=At 1,1) 

XMAX=AI N, 1 ) 

ABSCAI 1 )=XMIN 
ABSCAI 11)=XMAX 
OPDii'lI 1 )=YMIN 
URDINI 1 1) =YMAX 
DO 9 1=2,10 

2 = 1-1 

ABSCAI I )=I XMAX-XMIM)*Z/10.0+XMIN 

9 OPDINI : ) = IYKAX-Y!'.]N)=>''Z/10.0+YMIN 
PRINT 100, lORDIN'IJ) ,J = 1,11) 

STEPX= I XMAX-XiMlN ) /lOO .0 
KDELX=1 

KLINE=1 

LINE=1 

DC 26 IND=1,N 

IFINGRI D. EQ.O) GO TO 200 
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201 

10 

11 

12 

13 

14 

15 



16 

17 



kstep=line 

GO TO 201 

KbTEP=( A( INO, 1 J-XMIN) /STEPX+1.5 
OC 10 J=2,MP1 

TEMPY( J ) = A( IND, J )=«'100.0/( YMAX-YMIN)-YSHFT 

IF(KLINE-LINE) 12,12,18 

DO 13 1=2,100 

KAXiSd) = IDASH 

DO 14 1=1,101,10 

KAXIS( I ) = ISTAR 

IF(KSTEP-LINE) 15,15,17 

DC 16 1=2, MPl 

K = TEMPY( I ) + 1.5 

MM=IVAR ( I-l ) 

KAXIS(K) = CHAR(MM) 

PRINT 101, ABSCA( KDELX) , (KAXIS( J ) , J=l, 101 ) 
IF(MGRIO.EQ.O) GO TO 202 
KLINE=KLINE+NGRID 
ABSCA(KOELX) =A(KLINE, 1 ) 



GO TO 24 

202 KLINE=KLINE+10 
KDELX=K0ELX+1 
GO TO 24 

18 DC 19 1=2,100 

19 KAXIS( I ) = I BLANK 
DO 20 1=1,101,10 

20 KAXIS( I ) = IPER 
IF(KSTEP-LINE) 21,21,23 

21 DO 22 1=2, MPl 

K = TEMPY( I ) + 1.5 
MM=IVAk (1-1) 

22 KAXIS(K) = CHAR(MM) 

23 PRINT 102, (KAXIS( J),J=1,101) 

24 LINE=LINE+1 



IF (LINE-102) 25,25,27 

25 IF(KSTEP-LINE) 26,11,11 

26 CONTINUE 

27 RETURN 

end 
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